#!/usr/bin/python
"""
This program makes data points for GC skew plot

Usage:
Examples: 
dissertation_GCskew.py GKIL.v6.gbf percent
dissertation_GCskew.py GKIL.v6.gbf skew

Note: 

Author: Jimmy Saw
Date: 04-28-2012


"""

import sys
import re
import matplotlib.pyplot as plt
import pylab
import random
from Bio import SeqIO
from Bio.SeqUtils import GC

genome = SeqIO.read(sys.argv[1], "gb")
genome_size = len(genome.seq)

choice = str(sys.argv[2])

def calAvg(countlist):
    total = 0
    count = len(countlist)
    for i in countlist:
        total += i
    avg = float(total/count)
    return avg

def skew(seq):
    g = 0
    c = 0
    a = 0
    t = 0
    for i in seq:
        if i == 'G':
            g += 1
        elif i == 'C':
            c += 1
        elif i == 'A':
            a += 1
        elif i == 'T':
            t += 1
    gc = g + c
    at = a + t
    total = gc + at
    percentgc = (float(gc) / float(total)) * 100
    gcskew = float(g-c)/float(g+c) * 100
    #return gcskew
    if choice == 'skew':
        return gcskew
    elif choice == 'percent':
        return percentgc

#calculate skew
skewpoints = []

x = 1000
while x < len(genome.seq):
    start = x - 1000
    stop = x
    gc_skew = skew(genome.seq[start:stop])
    skewpoints.append((start, stop, gc_skew))
    x += 1000 #increment by 1kb

for i, j in enumerate(skewpoints):
    print "chr1", j[0], j[1], j[2]


